%% Main module of the projected power method for joint alignment proposed in the paper
% ``The Projected Power Method: An Efficient Algorithm for Joint Alignment from Pairwise Differences'' by Y. Chen and E. J. Candès.


function [zt, estErrors, z0, z0_time] = noncvxAlign(param, X_gt, L_opt, proj)
%  Projected power method for joint alignment
%    param:   parameters
%    X_gt:    ground truth (matrix representation)
%    L_opt:   multiplication operator ( L_opt(x) = L* x )
%    proj:    projection operator 

m = param.m;  n = param.n;  

%%  Spectral initialization
z0_time_start = tic;
% Compute rank-m approximation of L
U = randn(n*m, m);
B = L_opt( U );

% Block power method with the assistance of QR decomposition
for iter = 1: param.max_poweriter
    [U,Sigma]  = qr(B, 0);     % economy-size QR decomposition
    B = L_opt( U );
    
    if norm(B - U*Sigma) / norm(Sigma) <= param.eps_power, break;  end
end

% Round a random column to generate an initial guess
X0 = U * diag(diag(Sigma)) * U( randi(n*m), :)';
zt = proj( X0 ); 

z0 = zt;
z0_time = toc(z0_time_start);

% Error of initial guess
estErrors = computeError(zt, X_gt);
  

%% Projected power method
for iter = 2:  (param.maxIter - 1)
    zold = zt;
    zt   = L_opt( zt ); 
    zt   = proj( zt );
    estErrors = [estErrors, computeError(zt, X_gt)]; 
    
    if norm(zold - zt) == 0,  break;  end
end

% Round the final iterate to ensure feasibility

zt = ProjectOperator( zt, n, m, param.mu0 );
estErrors = [estErrors, computeError(zt, X_gt)]; 
 

%% Compute the L2 error of z modulo global offset
function dist = computeError(z, X_gt)

MSE = norm(z)^2 - 2* z'* X_gt + size(X_gt,1) / size(X_gt,2);  % || z - x ||^2 = || z ||^2 - 2 z'* x + n
dist = min( sqrt(MSE) );

